Monte Carlo simulations of the solid-liquid transition in hard spheres and 
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Monte Carlo simulations at constant pressure are performed to study coexistence and interfacial 
properties of the liquid-solid transition in hard spheres and in colloid-polymer mixtures. The latter 
system is described as a one-component Asakura-Oosawa (AO) model where the polymer's degrees 
of freedom are incorporated via an attractive part in the effective potential for the colloid-colloid 
interactions. For the considered AO model, the polymer reservoir packing fraction is rf v — 0.1 and 
the colloid-polymer size ratio is q = cr p /a = 0.15 (with <t p and a the diameter of polymers and 
colloids, respectively). Inhomogeneous solid-liquid systems are prepared by placing the solid fee 
phase in the middle of a rectangular simulation box creating two interfaces with the adjoined bulk 
liquid. By analyzing the growth of the crystalline region at various pressures and for different system 



sizes, the coexistence pressure p co is obtained, yielding p co = 11.576 



for the hard sphere system 



and Pco = 8.0 for the AO model (with ks the Boltzmann constant and T the temperature). 
Several order parameters are introduced to distinguish between solid and liquid phases and to 
describe the interfacial properties. From the capillary-wave broadening of the solid-liquid interface, 
the interfacial stiffness is obtained for the (100) crystalline plane, giving the values 7 « 0.49 ^f- 
for the hard-sphere system and 7 ta 0.95 for the AO model. 

PACS numbers: 



I. INTRODUCTION 

Various colloidal systems are ideal models for the in- 
vestigation of crystal nucleation and crystal growth pro- 
cesses. Whereas in atomistic systems, nucleation rates 
or interfacial free energies are hardly experimentally ac- 
cessible from direct measurements, in colloidal systems 
the much larger length and time scales allow to deter- 
mine these properties, at least in principle. For in- 
stance, in situ measurements of static structure factors 
in hard sphere-like colloidal systems using light scatter- 
ing techniques resulted in estimates of nucleation rates 
and have given insight into the applicability of classi- 
cal nucleation thcoryir— . Moreover, confocal microscopy 
gives a direct access to particle trajectories and thus, 
similar as in a computer simulation, any quantity of in- 
terest can be computed from the positions of the parti- 
cles. Recently, several experimental studies using con- 
focal microscopy^— have revealed various properties of 
solid-liquid interfaces. However, a direct estimate of 
anisotropic interfacial free energies in colloidal systems 
has not been possible so far. 

Two paradigms of colloidal model systems that 
can be realized experimentally are hard spheres and 
hard spheres with a short-ranged attractive interaction 
(colloid-polymer mixtures). Due to the short-range of 
the interactions, these model systems are also well- 
suited for theoretical studies (e.g. in the framework of 
density functional theory^) and for computer simula- 
tions. As far as the solid-fluid transition in hard spheres 
is concerned, various Molecular Dynamics (MD) and 
Monte Carlo (MC) techniques have been used to es- 
timate thermodynamic properties such as the coexis- 



tence pressure^—, kinetic growth coefficients^, and 
interfacial free energies^ - — . For the Asakura-Oosawa 
(AO) model of colloid-polymer mixtures (see below), MC 
studic o 19 ' 20 have provided estimates for the solid-liquid 
phase boundaries in a wide range of model parameters. 
However, to our knowledge, interfacial free energies for 
solid-fluid interfaces have not been determined so far for 
the AO model. 

Despite the efforts that have been recently undertaken, 
the examination of solid-liquid interfaces even for hard- 
sphere and hard-sphere-like systems is still subject to 
various problems. In particular, the role of finite-size 
effects has not been investigated in a systematic man- 
ner. Recently^, we have made a first preliminary step 
to fill this gap considering the solid-liquid interfaces of 
hard spheres and of the metallic system Ni. In the latter 
work, we have estimated the coexistence pressure and the 
interfacial stiffness 7 (see below) using constant-pressure 
MC simulations of solid-liquid inhomogeneous systems. 
To estimate 7, a result of capillary wave theory (CWT) 
was employed, according to which for a rough interface 
the mean-squared width of the interface grows logarith- 
mically with the lateral size of the system. Thus, we have 
made use of finite-size effects to compute the interfacial 
stiffness 7. 

In the present work, much more extensive calcula- 
tions are performed to determine the coexistence pres- 
sure and various interfacial properties. In addition 
to the hard-sphere system, the AO model for colloid- 
polymer mixtures is considered. As a matter of fact, it 
is much more difficult to compute interfacial properties 
for the latter colloid-polymer model than for the hard 
sphere system due to a slower growth kinetics as well as 
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smaller-amplitude capillary fluctuations along the inter- 
face. However, we show that both for the hard sphere 
system and the AO model accurate values for the coex- 
istence pressure are obtained without the requirement of 
extrapolation from relatively small system sizes to the 
thermodynamic limit. 

The rest of the paper is organized as follows. In the 
next section the interaction models are introduced and 
the main details of the simulation are given. The results 
on the solid-fluid coexistence pressure are worked out in 
Sec. IIII A[ followed by the introduction of local order pa- 
rameters and the discussion of the interfacial structure 
in Sec. IfflBl Then, Sec. iHTCl is devoted to the deter- 
mination of the interfacial stiffness from the finite-size 
broadening of the interface. Finally, a summary and dis- 
cussion of the results is provided in Sec. IIVI 
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II. MODEL SYSTEMS AND SIMULATION 
TECHNIQUES 

The one-component hard sphere model can be consid- 
ered as the simplest model with a solid-liquid transition. 
For a system of hard spheres of diameter er, the interac- 
tion potential is defined by 

with r the distance between two particles. The freezing 
of hard spheres has been first observed in early molecular 
dynamics simulations^—. In these simulations, systems 
of N = 500 particles were considered. In this work, sys- 
tems of up to 10 particles are investigated; in particular 
to obtain a reliable estimate of the interfacial stiffness 7. 

To a very good approximation, the equation of state 
for the fluid phase of the hard sphere model is given by 
the analytical Carnahan-Starling equatio n 24 ' 25 and to a 
lesser degree for the solid by the Hall equation^. Ac- 
curate values for the thermodynamic properties of the 
bulk solid are provided by empirical fits to computer 
simulations^. Since any allowed hard-sphere configu- 
ration has zero potential energy, the solid-fluid transi- 
tion in the hard-sphere system is completely driven by 
entropy and temperature T plays the role of a scaling 
factor. As a result, the thermodynamic properties are 
fully controlled by the packing density 77 = ^f-y (or, by 
the pressure p in case of fluctuating total volume V). At 
first glance, it seems to be surprising that hard spheres 
solidify since one may expect the entropy of an ordered 
solid phase to be always lower than that of the disordered 
fluid phase. However, at sufficiently high packing frac- 
tions, the spheres in a solid fee configuration have locally 
more freedom to move than in a fluid at the same pack- 
ing fraction, and the resulting higher number of possible 
microstates for the solid phase corresponds to a higher 
entropy. 

Experimentally, the repulsive interactions between 
hard sphere colloids can be modified by the addition of 



FIG. 1: A sketch of the phase diagram for the colloid-polymer 
mixture in the r/ P — r\ plane, assuming a polymer-colloid size 
ratio q = 0.15. Throughout this work, we consider either 
rjp — 0.0 (hard sphere model) or rf p — 0.1 (marked by the 
horizontal line). 



non-adsorbing polymers. Although the pairwise colloid- 
polymer as well as the polymer-polymer interactions are 
repulsive, an effective attraction between the colloids in 
colloid-polymer mixtures is induced entropically by a de- 
pletion effect^. 

A simple model for colloid-polymer mixtures is the 
Asakura-Oosawa (AO) model22r— where a hard sphere 
interaction is assumed between colloids as well as between 
colloids and polymers while polymer particles do not in- 
teract with each other. The limits of this approxima- 
tion are discussed elsewher o 33 ' 34 . In this two-component 
model, the strength of the attractive interactions between 
the colloids is controlled by the density of polymer parti- 
cles and the range of attraction by the size ratio between 
polymers and colloids, q = a p ja (with a p and a corre- 
sponding to the diameter of polymers and colloids, re- 
spectively). If one considers the system to be coupled to 
a polymer reservoir, the fugacity of polymers z p (or the 

polymer reservoir packing fraction r/ p = — |-^) can be re- 
garded as the analog of inverse temperature in a molecu- 
lar system and the phase diagram can be displayed in the 
?7p — 77 plane (corresponding to the temperature-density 
plane in a molecular system). At rf p = the coexistence 
region reduces to the case of pure hard spheres, whereas 
the increase of the polymer fugacity z p broadens the co- 
existence region, corresponding to a higher coexistence 
packing fraction of the solid phase and a lower one of the 
fluid phase (Fig.QJ. Note that throughout this work the 
polymer reservoir packing fraction is fixed at rf p = 0.1, 
as indicated by the horizontal line in Fig. [TJ 

On a qualitative level the AO model provides a good 
description of a colloid-polymer mixture. Modifications 
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of this model that lead to a more reliable description have 
been proposed elsewher e 33 ' 35 " —. On the other hand, the 
AO model can be simplified in terms of the computa- 
tional load by integrating out the polymers degrees of 
freedom and represent the colloid-polymer mixture as a 
one-component system of colloidal particles that inter- 
act with each other via an effective interaction potential. 
While for q > 0.154, this potential is a sum of two-body, 
three-body and higher-body terms, for q < 0.154 it re- 
duces to a pair potential given b y 30 ' 31 ' 38 
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a < r < a 
r > a + (T p 
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with r the distance between two colloids 
(fceT) -1 (with fee the Boltzmann constant). 

In this work, the fluid-to-solid transition of the hard 
sphere system and the AO model, as described by the 
potential $Z§, is studied using Monte Carlo (MC) sim- 
ulations in the isothermal-isobaric NpT and Np z T en- 
sembles, i.e. at constant pressure p, temperature T and 
particle number N (in the Np z T ensemble, p z is the di- 
agonal component of the pressure tensor perpendicular 
to the xy plane, i.e. perpendicular to the interface). MC 
simulations were carried out using a standard Metropolis 
algorithm. The trial moves were particle displacements 
and a rescaling of the volume for one MC cydeMr— . The 
maximum particle displacement was chosen such that the 
acceptance rate maintained constant at 30% for particle 
displacements and at 10% for volume's rescaling. To op- 
timize the speed of the simulation for large system sizes, 
a cell-linked neighbor list was use d 40 ' 41 . 



III. RESULTS 
A. Coexistence pressure 

As a prerequisite for the investigation of interfacial 
properties, an accurate determination of the coexistence 
pressure p co is required. To this end, inhomogeneous sys- 
tems are prepared where the crystal phase in the middle 
of an elongated simulation box is surrounded by the fluid 
phase, separated by two planar interfaces (note that two 
interfaces appear due to the use of periodic boundary 
conditions). When solid and liquid are in equilibrium at 
the pressure p co , the thermodynamic driving force is zero 
and the average total volume (V(t)) of the system does 
not change as a function of time, i.e. the crystal neither 
grows nor melts. In this manner, p co can be identified as 
the pressure where the time derivative (dV(t)/dt) van- 
ishes. 

To prepare an inhomogeneous system at a given pres- 
sure, one has to first compute separately the equation of 
state, p(rf), for the pure fluid and the pure crystal. From 
the fluid and the solid branch of p(rj) one can then read 
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FIG. 2: Phase diagram of the hard sphere system in the (p, rf) 
plane. The filled circles are simulation results, whereas the 
solid lines correspond to analytical expressions estimates of 
the equation of state, as proposed by Carnahan and Starling 
(fluid branch) and by Hall (solid branch). p co fts 11.576 is the 
estimated coexistence pressure at which the transition from a 
fluid to a crystalline fee phase occurs. The coexistence region 
is between the freezing point at rjt « 0.492 and the melting 
point at rjm ~ 0.545. 
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FIG. 3: Phase diagram of the AO model in the (p, rj) plane 
for q = 0.15 and rf p = 0.1. The filled circles are the simula- 
tion results. Here, the solid lines are spline interpolations be- 
tween simulation points. The estimated coexistence pressure 
IS Pco Rj 8.0. Freezing and melting points are at rft w 0.494 
and rj m w 0.64, respectively. 



off the packing fraction r/ (or the volume V) of the two 
phases at a given pressure. For the calculation of p(rj), we 
used MC simulations in the NpT ensemble. These runs 
were done for systems of N — 500 and, in some cases, 
for N = 1728 particles. The simulations with the larger 
system size indicated that finite-size effects are negligible 
for the calculation of the equation of state, at least for 
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FIG. 4: Static structure factors S(q) of the hard sphere and 
the AO fluid at the freezing point. 



systems with N > 500 particles. The simulation time 
was chosen depending on the convergence of the results, 
ranging from 500000 to several million MC cycles. 

Figures [2] and [3] display respectively the equation of 
state of the hard sphere and the AO systems (with 
q = 0.15 and rf p = 0.1). For the hard sphere system, 
the simulation results are well-described by analytical 
expressions (solid lines), as proposed by Carnahan and 
Starling^ for the liquid branch and by Hall2& for the fee 
phase. In the case of the AO system, no accurate ana- 
lytical predictions are available and so the solid lines in 
Fig. [3] are just spline interpolations that connect the data 
points. Also indicated in Figs. [5] and [3] are the coexistence 
pressure p co and the corresponding packing fractions at 
freezing and melting, rj = rjf and r\ = r] m , respectively, 
as obtained from our MC simulations (see below). For 
the AO model, the coexistence region is much broader 
than for the hard sphere system. However, the freezing 
point is at a similar packing fraction around r\ m 0.493 
for both models. Also the fluid structure at the freezing 
point is quite similar for the two systems. This can be 
inferred from Fig. @] where the static structure factor- 
s'^) is displayed for the hard sphere and the AO system 
at rj = rjf. From these findings different properties of the 
solid-fluid interface may be expected for the AO system 
when compared to the hard spheres. While at coexis- 
tence the fluid structure and fluid density are very sim- 
ilar for both model systems, in the AO model the fluid 
coexists with a fee crystal with a much higher packing 
density (rj m « 0.64) than in the hard sphere case where 
77m ~ 0.545. Below we see that indeed the interfacial 
stiffness for the AO system is about a factor of 2 higher 
than that of the hard spheres. 

Having computed the equations of state for the fluid 
and the solid phases, we now switch to the simulations 
of the inhomogeneous systems. As described above, the 
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FIG. 5: The relative change of colloidal packing fraction, Arj, 
as a function of Monte Carlo cycles for different pressures, 
as indicated. The solid lines are fits from which (dV/dt) is 
determined, a) Hard sphere system, b) AO model with q = 
0.15 and rf v — 0.1. In both cases, systems with iV = 10240 
particles are considered (corresponding to n = 8). 



coexistence pressure p co is obtained from the analysis of 
the change of the total volume as a function of time. 
The crystal grows or melts, dependent on the pressure 
at which the system is considered, and thus, the coexis- 
tence point is estimated as the state where the total vol- 
ume of the system does not change with time. Despite 
the simplicity of this analysis, there are several pitfalls 
when applying this method. First, interactions between 
the interfaces due to periodic boundary conditions have 
to be eliminated. Therefore, we have done test runs with 
various ratios of L z /L. As a result, L z = 5L was found 
to be an optimal choice for the avoidance of interaction 
effects between the interfaces. Second, finite-size effects 
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FIG. 6: Volume velocity (dV/dt) for different system sizes, as 
indicated, a) Hard sphere system, b) AO model. 



need to be quantified to obtain reliable estimates of var- 
ious properties at coexistence, such as the pressure or 
the interfacial stiffness. Therefore, we have considered 
various system sizes with N = 2500, 4320, 6860, 10240, 
20000, 67500, and 160000 particles. Third, the crystal- 
liquid interfaces have to be prepared such that no arti- 
ficial strains are generated in the crystalline region. In 
particular, the use of the isotropic NpT ensemble is not 
appropriate for the simulation of solid-liquid coexistence. 
The uniform volume moves in NpT simulations lead in- 
evitably to strains in the xy plane of the crystal, since the 
number of lattice planes cannot change in these moves. 
Thus, the lateral dimensions of the system need to be 
fixed such that they are commensurate with the chosen 
integer number of lattice planes (note that the lattice 
spacing at a given pressure is known from the bulk sim- 
ulations of the pure fee phase). Then, the volume is 
allowed to fluctuate in z direction (perpendicular to the 
interfaces), keeping the pressure in that direction (p z ) 
constant. 
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FIG. 7: Coexistence pressure as a function of system size 
(in terms of the inverse number of particles, 1/N) for the 
hard sphere system and the AO model. The dashed lines 
indicate the estimated coexistence pressures in the thermo- 
dynamic limit, p co = 11.576 for the hard sphere system and 
Pco = 8.0 for the AO model. For the hard sphere system, re- 
sults of other studies are included (from Errington^, Wilding 
and Braced, Poison and Frenketii, and Vega and Noya^). 



To generate crystal-fluid samples at various pressures, 
first independent solid and fluid samples were simulated. 
The solid was put into a simulation box with dimensions 
L x L x 3L, thereby aligning the (100) plane of the fee 
crystal perpendicular to the z-axis. The box length L 
was chosen such that it corresponds to the solid density 
at the considered pressure. At the same density, a start- 
ing configuration for the fluid was generated by putting 
the particles randomly in a box of size Lx Lx 2L. Then, 
the box dimensions of the fluid in x- and y-dircction were 
kept fixed and the fluid was equilibrated by a MC sim- 
ulation in the Np z T ensemble. After 10 6 MC cycles, 
solid and fluid samples were sufficiently thermalized and 
put together into a simulation box of approximate size 
L x L x 5L (applying periodic boundary conditions in all 
three spatial dimensions) , followed by further simulations 
in the Np z T ensemble. In the latter runs, the positions 
of the solid particles were fixed for the first 10 5 MC cycles 
to equilibrate the crystal-fluid interface without melting 
away the crystal due to an unfavorable local packing of 
particles in the interface region after matching the fluid 
slab with the solid slab. Then, the solid particles were 
released for the rest of the simulation. Finally, a set of 
short runs over 3 x 10 4 MC cycles were performed in a 
wide range of pressures to obtain at each pressure the 
total volume as a function of time. 

As an example, Fig. [S] shows the relative change of the 
colloidal packing fraction, At], at different pressures for 
systems of 10240 particles. At this system size, the lateral 



dimension is given by L — L x = L y 



na with n 



lattice planes in units of the lattice constant a (of course, 
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a changes as a function of pressure both for the AO and 
the HS system). In the following, we indicate the system 
size in terms of the number n, considering box geometries 
of nominal size L x L x 5L. We also note that the time 
t is measured in units of the number of Monte Carlo 
cycles; of course, it cannot be directly translated into 
a physical time, but, since we are not interested in the 
growth kinetics here, this does not matter in the present 
context. 

The slopes (dV/dt), averaged over 10 independent con- 
figurations, arc displayed in Fig. |6] for different system 
sizes, ranging from n = 5 to n = 20 for the hard sphere 
system and from n = 6 to n = 10 for the AO system. 
The values for p COl as estimated via the interpolation to 
(dV/dT) = for the different system sizes, are plotted 
in Fig. [7] as a function of the inverse number of particles 
1/N. Obviously, both for the HS and the AO model, 
finite-size effects arc small in the considered range of sys- 
tem sizes and we obtain p co = 11.576 ± 0.006^^ for 
the hard spheres and p co = 8.0 ± 0.026%£ for the AO 
model. The corresponding packing fractions for freezing 
and melting are respectively given by r/f = 0.492 and 
Tjm = 0.545 for the hard spheres and by r]{ = 0.494 and 
TJm = 0-64 for the AO model (see also Figs. [2] and [3]). 

Also included in Fig. [7] are estimates of p co for the 
hard sphere system, as obtained from other simulation 
studies^— . In these studies, the use of thermodynamic 
integration techniques as well as the phase switch MC 
method allowed only the consideration of relatively small 
system sizes and so these results lie below the dashed line 
in Fig. [7] that marks the estimate of p co in the thermo- 
dynamic limit, as obtained from our simulation. This 
indicates the advantage of the methodology used in this 
work: relatively large system sizes can be simulated and 
thus it is not necessary to perform extrapolations to the 
thermodynamic limit (at least not for the systems con- 
sidered here) that may easily lead to systematic errors in 
the estimate of the coexistence pressure. 
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FIG. 8: Order parameter distributions for the AO model 
(dashed lines) and the hard sphere system (solid lines) for 
the liquid and the fee phase, as indicated; a) §696, b) ty. 



with 



B. Local order parameters and interfacial structure 

Having determined p co for the hard sphere and the AO 
model, we can now investigate properties of solid-fluid in- 
terfaces in these systems at coexistence. First, order pa- 
rameters have to be identified that allow to distinguish 
between solid-like and fluid-like local order around a par- 
ticle. A class of order parameters that is well-suited for 
this purpose are local bond order parameters, Qi(i), as 
introduced by Steinhardt et alM. They are based on the 
expansion into spherical harmonics Yi m : 

( 1 \ 1/2 

Qi(i)= Uitt E l^-l 2 ( 3 ) 



Qlm{i) = i E YtmWij),^)) , (4) 

where ry is the distance vector between a pair of neigh- 
boring particles i and j, Zi is the number of neighbors 
within a given cut-off radius, and 0(rij) and (f>(rij) are the 
polar bond angles with respect to an arbitrary reference 
frame. 

A variant of the order parameters ((3]) has been put 
forward by ten Wolde et by introducing the dot 

product 



QlQl(i) = ^"E E Qlrn(i)qim(j)* > ( 5 ) 
1 j — 1 m— — l 
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b) 9 for the AO model. The system contains 20000 particles, 
corresponding to n = 10. 
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Here, we use qeqe that is defined by Eqs. ([5]) and ([6]), 
setting I — 6. 

Figure|5£i shows the q^qe distributions for the pure fluid 
and fee phases at coexistence for the hard sphere system 
(solid lines) and the AO mixture (dashed lines). The 
relatively small overlap of the distribution for the solid 
with that of the fluid indicates that q§q§ is well-suited to 
distinguish between solid and fluid particles. Note that 
we have used time-averaged particle positions for the cal- 
culation of the order parameter distributions in Fig. [8^ 
(also for the distributions shown in Fig.|5b). Particle po- 
sitions were averaged over 50 MC cycles in case of the 
hard sphere system and over 20 MC cycles in case of the 
AO mixture. This reduces the shift of the order param- 
eter distributions for the solid fee phases to lower values 



of the order parameter, as compared to the distribution 
of an ideal fee crystal. 

A different local order parameter has been introduced 
by Morri a 44 ' 45 . It is defined by 



1 1 Z i N 1 

—-YY 

q 3=1 fc=l 



cxp(iq k ■ fij) 



(7) 



where denotes the distance vector of a particle i to 
neighboring particles j, and the wave- vectors q% are re- 
lated to reciprocal vectors of the fee lattice with lattice 
constant ao, <fi = 27r/ao(— 1 , 1, —I), q\ = 27r/arj(l, — 1, 1) 
and q\ — 27r/ao(l, 1,-1). An additional average of ^(i) 
over a particle with index i and its neighboring particles 
in the first coordination shell yields 



*(*) 




(8) 



As can be inferred from Fig. [SJ}, the ^(i) distributions 
display a sharp peak close to $ = 1 for the fee phase and 
one close to *f> = for the fluid. Obviously, the order 
parameter ^> is also well-suited to identify the local order 
around particles. 

By applying the methodology of our previous study^, 
we characterize the local structure of the interfaces by z- 
dependent profiles of averaged local order parameters. 
For this purpose solid-fluid samples were divided into 
bins of length Az = 0.05c along the z-direction and 
"time" averages of the order parameter were obtained 
for each bin, thereby correcting for shifts of the crys- 
tal planes with respect to the reference frame during the 
simulation. 

Figures [5Ji and [SJd show order parameter profiles for 
the AO model as a function of the z-coordinate (i.e. per- 
pendicular to the interface). To obtain these profiles, 
the sum of the order parameter of the particles in each 
bin transversal to the solid-fluid interface is divided by 
the volume of the bin AV^ = L 2 Az. In the crystalline 
region, the profiles exhibit strong oscillations with the 
periodicity of the crystalline planes. The amplitude of 
these oscillations decay rapidly in the interfacial region 
and flatten completely in the fluid region. Correspond- 
ing data for the hard-sphere system are reported in our 
previous studyi^. 



Finite-size interfacial broadening and interfacial 
stiffness 



Now, we aim at computing the interfacial stiffness 7 
which is defined by 7 = 7 + d 2 7 /d8 2 with 7 the interfa- 
cial tension and 9 the angle between the interface normal 
and the (100) direction. Whereas 7 describes the free 
energy cost to increase the area of the interface, the term 
cP^/dO 2 accounts for the free energy required to locally 
change the orientation of the crystal. The latter term 
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FIG. 10: Coarse-grained order parameter profiles for the hard 
sphere system for different system sizes; a) q^qe, b) ty. The 
lines are fits to Eq. ([9]). 



FIG. 11: Coarse-grained order parameter profiles for the AO 
model for different system sizes; a) qaqe, b) The lines are 
fits to Eq. ©. 



would of course vanish for a system where both phases 
at the interface are isotropic (as, e.g., in the case of liquid- 
gas interfaces). 

In the framework of capillary wave theory (CWT), the 
interfacial stiffness can be computed from the broaden- 
ing of the solid-fluid interface with increasing lateral size 
of the system^ (see below). To make use of this predic- 
tion, the apparent mean-squared width of the interface, 
w 2 , has to be determined and analyzed as a function of 
system size. For this purpose, fine-grid order profiles as 
the ones shown in Fig. [S] are not well-suited since the os- 
cillations due to the crystalline order do not allow to fit 
the profile with a simple function where the interfacial 
width appears as a free parameter. Therefore, it is useful 
to coarse-grain the profiles by averaging over the oscilla- 
tions. To this end, we have identified the minima in the 
fine-grid profiles and used them to mark the borders of 
non-uniform bins that match the crystalline layers. Then, 
we took the average of the order parameter in each of the 



latter bins. 

Examples of the resulting coarse-grained profile for dif- 
ferent system sizes are shown in Figs. QJJ] and [IT] As we 
can see in these plots, the data can be well fitted with 
hyperbolic tangent function of the following form 

,. , A + B A-B /z-Zo\ 
*( z ) = _ + _tanh^— -J (9) 

where A and B are the bulk order parameters in the 
solid and fluid obtained from the independent bulk sim- 
ulations. So the position of the interface zq and its width 
w are the only fit parameters. Note that Eq. (|9]) can be 
obtained in the framework of mean-field theory. How- 
ever, here it is just used a fitting function to determine 
the width w as a function of the lateral size of the sys- 
tem. Figures [10] and [TT] indicate that both for the hard 
sphere model and the AO mixture, w slightly increases 
with increasing system size. 

The finding that the width w is system-size-dcpcndcnt 
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a) 



HS 




b) 



AO 




FIG. 12: Snapshots of the interfaces for a) hard spheres and 
b) the AO model for systems with n = 15. Here, only solid 
particles are shown, defining a particle as a solid one when its 
qaqa > 0.6 and its coordination number greater than 10. 



is expected from CWT where the mean-squared interfa- 
cial width is found to diverge logarithmically with the 
lateral size L of the system^ - —, 



w 



w 2 



k B T. L 
In— 

4 7 I 



(10) 



with wq the so-called intrinsic width and / a cut-off 
parameter that has to be introduced since CWT only 
describes long-wavelength undulations of the interface. 
Note that it is impossible to disentangle the intrin- 
sic width contribution Wg from the cut-off contribution 
■^j^lnZ. An extensive discussion of this issue has been 

provided for the case of polymer mixtures^ - — . 

The existence of capillary waves is associated with 
the roughness of the interface and the applicability of 
CWT requires that interfaces above a possible rough- 
ening transition^— are considered. That both for the 
hard sphere system and the AO mixture the solid-fluid 
interface is indeed rough, is indicated by the snapshots in 
Fig. [HI Here, only the solid particles are shown, defining 
a particle as a solid one if qeq& > 0.6 and if its nearest- 
neighbor coordination number is larger than 10 (note 
that we have only needed this definition for the snapshots 
in Fig. IT2|) . The snapshots reveal that the solid-fluid in- 
terface of the hard sphere system is significantly "more 
rough" than that of the AO mixture. This might be due 
to the larger density difference between the fee and the 
fluid phase in case of the AO model, when compared to 
the hard sphere system. 

In Fig. [T3J the prediction (JTUJ) is confirmed: both for 
the hard spheres and the AO mixture, the mean-squared 
width w 2 exhibits a logarithmic divergence with respect 
to the lateral dimension L. Within the statistical errors, 
the two order parameters q^q^ and <F yield similar results. 
From the fits with Eq. (fTU|) . we estimate for the (100) 
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FIG. 13: Mean-squared width w 2 as a function of InL, for 
a) hard spheres and b) the AO model. The values for 7 are 
obtained from the fits (solid lines) to Eq. (|10|l . 



orientation 7 w 0.49 ± 0.02kBT/cr 2 for the hard sphere 
system and 7 w 0.95 ± Q.WbT/o 2 for the AO mixture. 
The value of 7 for the hard sphere system roughly agrees 
with previous estimates, obtained by other methods (for 
a discussion of this issue, see Refsi 18 ' 60 ). 



IV. CONCLUSIONS 

We have investigated the fluid-to-crystal transition of 
two paradigmatic colloidal model systems using Monte 
Carlo (MC) simulations at constant pressure. Inhomo- 
geneous systems have been prepared where the crystal 
phase in the middle of an elongated simulation box is 
separated from the fluid phase by two interfaces. We 
have demonstrated that MC simulations with this setup 
allow for reliable estimates of the coexistence pressure, 
p co , and the interfacial stiffness, 7. Both for the hard 
sphere system and the AO mixture, our methodology al- 
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lows for the simulation of relatively large systems and 
thus, the coexistence pressure p co can be computed with- 
out relying on error-prone extrapolations to the thermo- 
dynamic limit. On the other hand, we have presented a 
method for the calculation of 7 that makes use of finite- 
size effects due to the capillary-wave broadening of the 
interface. So far, this method has been mainly applied to 
polymer interfaces^, liquid-vapor interfaces^, isotropic- 
nematic interfaces^ or interfaces in Ising systems^, but 
not to solid-fluid interfaces, as in this work (an exception 
is of course our recent preliminary work on hard spheres 
and nickel^). It requires a systematic variation of the 
system size and relatively small statistical errors to re- 
solve the logarithmic divergence of the mean-squared in- 
terfacial width, as predicted in the framework of capillary 
wave theory (CWT). With respect to the latter issues, 
the determination of 7 is much more difficult for the AO 
mixture than for the hard sphere system. First, in the rel- 
evant range of colloid packing fractions rj mass transport 
processes of the AO fluid in the bulk and in the interface 
region are much slower than in the hard sphere system. 
Therefore, much longer simulation runs are necessary for 
the AO model, to achieve comparable statistics as in the 
hard sphere system. Second, the interfacial stiffness 7 
for the AO model is about a factor of two higher than in 
the hard sphere system, associated with lower-amplitude 
capillary fluctuations at comparable system sizes. Thus, 
the signal-to-noise ratio is worse for the AO model when 
one analyzes the interfacial fluctuations. However, de- 
spite the latter difficulties, the estimate of 7 from the 



interfacial broadening works also well for the AO model. 

CWT can be of course only applied to the analysis of 
interfacial properties, if rough, non-faceted interfaces are 
considered. Otherwise, there would be no divergence of 
the width of the interface with increasing lateral size of 
the system. For the case of the model systems consid- 
ered in this work, the applicability of CWT is justified 
since we observe the logarithmic growth of the mean- 
squared interfacial width as a function of the lateral di- 
mension L. To further rationalize the use of CWT, we 
plan to study the Fourier spectrum of interfacial fluctu- 
ations from which one can alternatively determine the 
interfacial stiffness. Moreover, we plan to investigate the 
possibility of a roughening transition of the AO model 
where the solid-fluid interface changes from a faceted to a 
rough interface with long- wavelength capillary wave fluc- 
tuations. 
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